function hessianEval()
% Get hessian matrix of NLP
%
global prob nlp
global hstruct

% global sigma lambda
% sigma = 1; lambda = ones(nlp.ncon, 1);
% if strcmp(nlp.hessian, 'autodiff') == 1
%     hpi = gethpi('lagAD', nlp.nvar);
%     hstruct = tril(hpi.SPH);
%     % Set function handle
%     nlp.func.hessian = @hessianAD;
%     nlp.func.hessianstructure = @() hstruct;
% end

global iG iG_cut iGB_cut iLD Bgen_cut Gld_cut Bld_cut
global EXd EXdTJ
global rows cols aux
global iH jH vH

if strcmp(nlp.hessian, 'analytic') == 1
    calcBaseHessian();
    calcAuxiliary();
    nlp.func.hessian = @hessianAN;
    nlp.func.hessianstructure = @() hstruct;
end

%% sub-functions
% 	function H = hessianAD (x, sgm, lmbd)
% 		sigma = sgm; lambda = lmbd;
% 		[~,~,H] = evalh('lagAD', x, [], hpi);
%         H = tril(H);
%     end
    function H = hessianAN (x, ~, lmbd)
        nnode = nlp.nnode; nG = prob.nG; 
        ncutG = prob.ncutG; ncutL = prob.ncutL;
        uGa = repmat(x(nlp.uGidx),1,nnode)'; 
        EXd1uG = aux.EXdc .* (1-uGa(:));
        lB = lmbd(aux.lidxB); 
        lCgx = lmbd(aux.lidxCgx); lCgy = lmbd(aux.lidxCgy); 
        lClx = lmbd(aux.lidxClx); lCly = lmbd(aux.lidxCly); 
        dlt = x(aux.hdltidx); dltc = x(aux.hdltcidx);
        Vx = x(aux.hVxidx); Vy = x(aux.hVyidx);
        Vxg = x(aux.hVxgidx); Vyg = x(aux.hVygidx);
        % <1> D_dlt_dlt
        cache1 = aux.EXdTJ.*(-Vxg.*sin(dlt)+Vyg.*cos(dlt)) .* lB;
        cache11 = EXd1uG.*sin(dltc) .* lCgx ...
            - EXd1uG.*cos(dltc) .* lCgy; % EXd1uG: EXd.*(1-uG)
        cache1(aux.hdltcloc) = cache1(aux.hdltcloc) + cache11;
        vidx = 1 : nnode * nG; vH(vidx) = cache1;
        % <2> D_Vx_dlt
        cache1 = aux.EXdTJ.*cos(dlt) .* lB;
        vidx = vidx + nnode * nG; vH(vidx) = cache1;
        % <3> D_Vy_dlt
        cache1 = aux.EXdTJ.*sin(dlt) .* lB;
        vidx = vidx + nnode * nG; vH(vidx) = cache1;
        % <4> D_uG_dltc
        cache2 = aux.EXdc.*cos(dltc) .* lCgx + aux.EXdc.*sin(dltc) .* lCgy;
        vidx = 3*nnode*nG +(1:nnode*ncutG); vH(vidx) = cache2;
        % <5> D_uG_Vcx
        vidx = vidx + nnode*ncutG; vH(vidx) = -aux.Bg .* lCgy;
        % <6> D_uG_Vcy
        vidx = vidx + nnode*ncutG; vH(vidx) = aux.Bg .* lCgx;
        % <7> D_uL_Vcx
        vidx = 3*nnode*(nG + ncutG) +(1:nnode*ncutL);
        vH(vidx) = -aux.Gl .* lClx - aux.Bl .* lCly;
        % <8> D_uL_Vcy
        vidx = vidx + nnode*ncutL;
        vH(vidx) = aux.Bl .* lClx - aux.Gl .* lCly;
        % build Hessian
        H = sparse(iH, jH, vH, nlp.nvar, nlp.nvar);
    end
    % [iH jH] of Hessian matrix
    function calcBaseHessian()
        nnode = nlp.nnode; nstate = prob.nstate; 
        nG = prob.nG; nB = prob.nB; ncutG = prob.ncutG; ncutL = prob.ncutL;
        nzH = 3*nnode*(nG + ncutG) + 2*nnode*ncutL;
        iH = zeros(nzH, 1); jH = zeros(nzH, 1); vH = zeros(nzH, 1);
        % <1-BC> D_dlt_dlt [NZ = nnode*nG] (3 types)
        idxH = 1 : nnode*nG; tmp = cols(:, 1:nG); cs = tmp(:);
        iH(idxH) = cs; jH(idxH) = cs;
        % <2-B> D_Vx_dlt  [NZ = nnode*nG] (1 type)
        idxH = idxH + nnode*nG; tmp = cols(:, nstate +iG); rs = tmp(:);
        iH(idxH) = rs; jH(idxH) = cs;
        % <3-B> D_Vy_dlt  [NZ = nnode*nG] (1 type)
        idxH = idxH + nnode*nG; rs = rs + nB*nnode;
        iH(idxH) = rs; jH(idxH) = cs;
        % <4-C> D_uG_dltc [NZ = nnode*ncutG] (2 types)
        e = idxH(end); idxH = e+1 : e+nnode*ncutG; 
        tmp = cols(:, iG_cut); cs = tmp(:);
        tmp = repmat(nlp.uGidx, nnode, 1); rs = tmp(:);
        iH(idxH) = rs; jH(idxH) = cs;
        % <5-C> D_uG_Vcx  [NZ = nnode*ncutG] (1 type)<const>
        idxH = idxH + nnode*ncutG; tmp = cols(:, nstate +iGB_cut); cs = tmp(:);
        iH(idxH) = rs; jH(idxH) = cs;
        % <6-C> D_uG_Vcy  [NZ = nnode*ncutG] (1 type)<const>
        idxH = idxH + nnode*ncutG; cs = cs + nB*nnode;
        iH(idxH) = rs; jH(idxH) = cs;
        % <7-C> D_uL_Vcx  [NZ = nnode*ncutL] (2 types)
        e = idxH(end); idxH = e+1 : e+nnode*ncutL;
        tmp = repmat(nlp.uLidx, nnode, 1); rs = tmp(:);
        tmp = cols(:, nstate +iLD); cs = tmp(:);
        iH(idxH) = rs; jH(idxH) = cs;
        % <8-C> D_uL_Vcx  [NZ = nnode*ncutL] (2 types)
        idxH = idxH +nnode*ncutL; cs = cs + nB*nnode;
        iH(idxH) = rs; jH(idxH) = cs;
        % build structure
        hstruct = sparse(iH, jH, 1, nlp.nvar, nlp.nvar);
    end
    % indices of Hessian
    function calcAuxiliary ()
        nnode = nlp.nnode; nstate = prob.nstate; 
        nG = prob.nG; nB = prob.nB;
        % lambda index
        aux.lidxB = nG*nnode+1 : nstate*nnode;
        tmp = rows(:,2*nG+iGB_cut); aux.lidxCgx = tmp(:);
        aux.lidxCgy = aux.lidxCgx + nB*nnode;
        tmp = rows(:,2*nG+iLD); aux.lidxClx = tmp(:);
        aux.lidxCly = aux.lidxClx + nB*nnode;
        % var index
        tmp = cols(:,1:nG); aux.hdltidx = tmp(:);
        tmp = cols(:,iG_cut); aux.hdltcidx = tmp(:);
        tmp = rows(:,iG_cut); aux.hdltcloc = tmp(:);
        tmp = cols(:, nstate+iGB_cut); aux.hVxidx = tmp(:);
        tmp = cols(:, (nstate+nB)+iGB_cut); aux.hVyidx = tmp(:);
        tmp = cols(:, nstate+iG); aux.hVxgidx = tmp(:);
        tmp = cols(:, (nstate+nB)+iG); aux.hVygidx = tmp(:);
        % expanded var
        cache0 = repmat(EXd, 1, nnode)'; aux.EXd = cache0(:);
        cache0 = repmat(EXd(iG_cut), 1, nnode)'; aux.EXdc = cache0(:);
        cache0 = repmat(EXdTJ, 1, nnode)'; aux.EXdTJ = cache0(:);
        cache0 = repmat(Bgen_cut, 1, nnode)'; aux.Bg = cache0(:);
        cache0 = repmat(Bld_cut, 1, nnode)'; aux.Bl = cache0(:);
        cache0 = repmat(Gld_cut, 1, nnode)'; aux.Gl = cache0(:);
    end
end